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Many problems at the forefront of theoretical astrophysics require the treatment of magnetized 
fluids in dynamical, strongly curved spacetimes. Such problems include the origin of gamma-ray 
bursts, magnetic braking of differential rotation in nascent neutron stars arising from stellar core 
collapse or binary neutron star merger, the formation of jets and magnetized disks around newborn 
black holes, etc. To model these phenomena, all of which involve both general relativity (GR) 
and magnetohydrodynamics (MHD), we have developed a GRMHD code capable of evolving MHD 
fluids in dynamical spacetimes. Our code solves the Einstein-Maxwell-MHD system of coupled 
equations in axisymmetry and in full 3-1-1 dimensions. We evolve the metric by integrating the 
BSSN equations, and use a conservative, shock-capturing scheme to evolve the MHD equations. 
Our code gives accurate results in standard MHD code-test problems, including magnetized shocks 
and magnetized Bondi flow. To test our code’s ability to evolve the MHD equations in a dynamical 
spacetime, we study the perturbations of a homogeneous, magnetized fluid excited by a gravitational 
plane wave, and we hnd good agreement between the analytic and numerical solutions. 
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I. INTRODUCTION 

Magnetic fields play a crucial role in determining the 
evolution of many relativistic objects. In any highly con¬ 
ducting astrophysical plasma, a frozen-in magnetic field 
can be amplified appreciably by gas compression or shear. 
Even when an initial seed field is weak, the field can grow 
in the course of time to significantly influence the gas dy¬ 
namical behavior of the system. In problems where the 
self-gravity of the gas can be ignored, simulations can be 
performed without numerically evolving the spacetime 
metric. Some accretion problems fall into this category. 
In many other problems, the effect of the magnetized 
fluid on the metric cannot be ignored, and the two must 
be evolved self-consistently. The final fate of many of 
these relativistic astrophysical systems, and their distin¬ 
guishing observational signatures, may hinge on the role 
that magnetic fields play during the evolution. Some of 
these systems are promising sources of gravitational radi¬ 
ation for detection by laser interferometers such as LIGO, 
VIRGO, TAMA, GEO and LISA. Some also may be re¬ 
sponsible for gamma-ray bursts. Examples of astrophysi¬ 
cal scenarios involving strong-held dynamical spacetimes 
in which MHD effects may play a decisive role include 
the following: 

• The merger of binary neutron stars. The coales¬ 
cence can lead to the formation of a Iwpermassive star 
supported by differential rotation P, Pi- While such a 
star may be dynamically stable against gravitational col¬ 
lapse and bar formation, the radial stabilization due to 
differential rotation is likely to be temporary. Magnetic 
braking and viscosity combine to drive the star to uni- 
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form rotation, even if the seed magnetic field and the 
viscosity are small P. This process can lead to delayed 
collapse and massive disk formation P , accompanied by 
a delayed gravitational wave burst. MHD-related insta¬ 
bilities in differentially rotating stars might also drive a 
gamma-ray burst P. 

• Core collapse in a supernova. Core collapse may 
induce differential rotation, even if the rotation of the 
progenitor at the onset of collapse is only moderately 
rapid and almost uniform (see, e.g. P, and references 
therein). Differential rotation can wind up a frozen-in 
magnetic Held to high values, at which point it may pro¬ 
vide a signihcant source of stress, which could affect the 
explosion P. 

• The generation of gamma-ray bursts (GRBs). Short- 

duration GRBs are thought to result from binary neutron 
star mergers Ip, or tidal disruptions of neutron stars by 
black holes P, or hypergiant Hares of ‘magnetars’ as¬ 
sociated with the soft gamma-ray repeaters Long- 
duration GRBs likely result from the collapse of rota ting , 
massive stars which form black holes (‘collapsars’) [Tl|. 
In current scenarios, the burst is powered by the extrac¬ 
tion of rotational energy from the neutron star or black 
hole, or from the remnant disk material formed around 
the black hole 0. Strong magnetic Helds provide a 
likely mechanism for extracting this energy on the re¬ 
quired timescale and driving collimated GRB outHows 
in the form of relativistic jets Even if the initial 

magnetic Helds are weak, they can be ampHHed to the 
required values by differential motions HI- 

• Supermassive black hole (SMBH) formation. The 
origin of the SMBHs observed in galaxies and quasars is 
one of the great mysteries of contemporary astrophysics. 
Several hypotheses for the origin of SMBHs involve rel¬ 
ativistic, self-gravitating fluids in which magnetic fields 
can play an important role. It is thought that SMBHs 
start from smaller initial seed black holes, which grow 
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to supermassive size by a combination of accretion and 
mergers. The seed black holes might be provided by 
the collapse of massive (M ~ 10 ^Mq) Population III 
stars If so, magnetic forces will affect the collapse 
leading to the formation of these seeds, as well as their 
growth by accretion [3113 Another possibility is that 
SMBHs form directly from the catastrophic collapse of 
supermassive stars (SMSs) 0 This collapse will pro¬ 
ceed differently, depending on whether the SMS rotates 
uniformly or differentially |l 9|| . Magnetic fields and tur¬ 
bulence provide the principle mechanisms that can damp 
differential rotation in such stars [^ and thus determine 
their ultimate fate. 

• The r-mode instability in rotating neutron stars. 
This instability has been proposed as a possible mech¬ 
anism for limiting the angular velocities in neutron stars 
and producing observable quasi-periodic gravitational 
waves [ 23 . However, preliminary calculations suggest 
that if the stellar magnetic field is strong enough, r-mode 
oscillations will not occur . Even if the initial field is 
weak, fluid motions produced by these oscillations may 
amplify the magnetic field and eventually distort or sup¬ 
press the r-modes altogether. (R-modes may also be sup¬ 
pressed b y no n-linear mode coupling or hyperon bulk 
viscosity j24l|.l 

• Massive disk accretion. The importance of magnetic 
effects on gas accretion onto a black hole is well known. In 
many cases, the density of the accreting material is small 
enough that its effect on the spacetime geometry is neg¬ 
ligible. Such systems can be studied by evolving the gas 
on the stationary Kerr spacetime background produced 
by the central black hole. There are, however, situa¬ 
tions in which accretion disks with masses comparable to 
that of the central black hole can be formed. Examples 
include the collapse of rapidly rotating stars or super¬ 
massive stars 1^, and neutron star merger (especially 
when the two neutron stars have unequal masses |3). In 
these cases, the spacetime is dynamical and Einstein’s 
equations for the metric must be evolved along with the 
MHD equations. 

In the recent years, numerical codes which evolve 
the general relativistic MHD equations on fixed 
Schwarzschild or Kerr black hole spacetimes have been 
developed by Yokosawa |2^, Koide et al [^, Komis¬ 
sarov [^, De Villiers and Hawley j^, and Gammie et 
al [31 • These codes have been used to stud y th e struc¬ 
ture of accretion flows onto Kerr black holes 13l j the 
Blandford-Znajek effect in low-density regions near the 
hole |3IS3i and the formation of GRB jets [3|- 

In contrast to the above effort, few attempts have been 
made to simulate relativistic MHD flows in dynamical 
spacetimes. One major attempt was by Wilson, thirty 
years ago |3l- He simulated the collapse of a rotating 
SMS with a frozen-in poloidal magnetic field in axisym- 
metry using a code which assumed the conformal flatness 
approximation for the spatial metric. No gravitational 
radiation is allowed in this approximation. Wilson’s work 
was generalized by Nakamura, Oohara and Kojima in 


1987 [3|- They studied the effect of poloidal magnetic 
fields on the collapse of nonrotating SMSs in full GR in 
axisymmetry. Since the stars are nonrotating, toroidal 
fields are not generated, which simplifies the calculation. 
Since then, no other simulations of this kind have been 
attempted. However, in anticipation of future numeri¬ 
cal work, formulations of the coupled Einstein-Maxwell- 
MHD equations were proposed by Sloan and Smarr [^ . 
by Zhang |3| , and by Baumgarte and Shapiro [3 • 

In this paper, we present the first code capable of evolv¬ 
ing the Einstein-Maxwell-MHD equations without ap¬ 
proximation in both two dimensions (axisymmetry) and 
three dimensions. Our code is based on the Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN) formulation of the 
3-1-1 Einstein field equations |4r| . In previous papers, 
we have evolved the BSSN equations coupled to a per¬ 
fect fluid ^3> have applied our code to simulate 

stellar collapse and binary neutron star inspiral ^3 ■ 
We then gmeralized our code to study fluids with shear 
viscosity Q, and we implemented black hole excision 
techniques to study the collapse of fluid stars to black 
holes |3l- In this paper, we have completely reformu¬ 
lated the hydrodynamics sector of our code in order to 
improve its accuracy and shock-handling capability. We 
have extended the code to allow for a magnetic field 
frozen into a perfectly conducting fluid in the MHD ap¬ 
proximation. 

In Sec. H, we present the evolution equations inte¬ 
grated by our code. In Sec. HI, we discuss the adopted 
numerical techniques. We perform a variety of code tests 
in Sec. IV, including magnetized shocks, a magnetized 
Bondi accretion flow, and a linear gravitational wave that 
excites MHD waves. 


II. FORMALISM 

A. Evolution of the gravitational fields 

Throughout this paper, Latin indices denote spatial 
components (1-3) and Greek indices denote spacetime 
components (0-3). We write the metric in the form 

ds^ = —a^dt^ + ')ij{dx'' + f3''dt){dx^ + fd^dt), (1) 

where a, 13'^, and "fij are the lapse, shift, and spatial met¬ 
ric, respectively. The extrinsic curvature Kij is defined 
by 

{dt - = -2aKij, ( 2 ) 

where is the Lie derivative with respect to /3L We 
adopt geometrized units, so that G = c = 1. We evolve 
7 ij and Kij using the BSSN formulation ^3- The fun- 
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damental variables for BSSN evolution are 


(t> 

= 

^ln[det(7*j)] , 

(3) 

% 

= 

e '^^lij , 


(4) 

K 

= 



(5) 

Aij 

= 


^lijE) , 

(6) 

r 

= 

-f\j ■ 


(7) 


The evolution and constraint equations for these fields 
are summarized in |4ft l43| . In the presence of mass- 
energy, these evolution equations contain the following 
source terms: 

P = nanf3T°‘^ , 

Si = , ( 8 ) 

Sij = , 

where T“^ is the stress tensor, and n“ = 
is the time-like unit vector normal to the t = constant 
time slices. Note that, since the electromagnetic fields 
contribute to T“^, they will contribute to p, Si, and Sij 
as shown in Eq. 

In order to evolve the 3-1-1 Einstein equations forward 
in time, one must choose lapse a and shift /3* functions, 
which specify how the spacetime is foliated. The lapse 
and shift must be chosen in such a way that the total 
system of evolution equations is stable. It is also desirable 
that the adopted gauge conditions make the evolution 
appear as stationary as possible. As in use the 

following hyperbolic driver conditions: 

dta = aA 

dtA = —ai{adtK ( 9 ) 

+a2dta + a-ie~'^‘^aK) . 

= b^{adtr-h2dtP^), ( 10 ) 

where ai, 02 , 03 , bi, and 62 are freely specifiable con¬ 
stants. We usually choose ai = 61 = 0.75, 02 = 62 = 
0.34M“^, 03 = 1. (There are exceptions. Eor example, 
when evolving a collapsing system, it is better to use a 
smaller bi. This prevents “blowi^ out” of the coordi¬ 
nate system, a well-known effect which can spoil 

grid resolution in the center of the collapsing object.) 

After the initial time, we do not enforce the constraint 
equations, but rather monitor them as a check on the ac¬ 
curacy of our evolution of the metric. Another check on 
the metric evolution is the conservation of the ADM mass 
M and angular momentum J of the spacetime, account¬ 
ing for losses due to gravitational radiation or rest-mass 
outflow from the grid (both are negligible for the tests re¬ 
ported here). The formulae for M and J are given in 
in terms of volume integrals over the whole space. Eor 
spacetimes containing black holes in which excision is em¬ 
ployed, we use the formulae given in |43|. which consist 
of a surface integral over a sphere enclosing the excision 
region, plus a volume integral from the surface to spatial 
infinity. 


B. Evolution of the electromagnetic fields 

The electromagnetic stress-energy tensor is given 

by 

TZ = ^ . ( 11 ) 

We decompose the Earaday tensor F^^" as 

F^^'^ = - rfE^ + , (12) 

so that E^^ and are the electric and magnetic fields 
measured by a normal observer n^. Both fields are purely 
spatial {E^Ufj, = B^n^ = 0), and one can easily show 
that 

E^^ = , (13) 

where 

F*^’' = (14) 

is the dual of F^^. In terms of E^ and B^, the electro¬ 
magnetic stress tensor is given by 

TZ = + 2nf^n‘'){E^Eg,+ B^Bx) 

Stt 

—L + E^^E'') 

47r 

^LnA^^'A^E^Bs . (15) 

47r 

Along with the electromagnetic field, we also assume the 
presence of a perfect fluid with rest density po? pressure 
P, and 4-velecity u^, so that the total stress-energy ten¬ 
sor is 

= pohu^u’' + Pg>^'' + TZ , (16) 

where the specific enthalpy h is related to the specific 
internal energy e by /i = 1 -I- e -I- F/po- 

Eor most applications of interest in relativistic astro¬ 
physics, we can assume perfect conductivity. In the limit 
of infinite conductivity, Ohm’s law yields the MHD con¬ 
dition: 

= 0 . (17) 

The electric and magnetic fields measured by an observer 
comoving with the fluid are [cf. Eq. (ICT l 

E^^^ = , B[^^ = . (18) 

The ideal MHD condition 03 is equivalent to the state¬ 
ment that the electric field observed in the fluid’s rest 
frame vanishes = 0). Note that is orthogonal 
to Ufj_, i.e. = 0. We can express F>^'' in terms of 

-^(1) 8'S [cf- Ecf- 03] 

F^^’' = . (19) 
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Taking the dual of Eq. (P|l. we obtain 

. ( 20 ) 

We define the projection operator = g^i, + u^Uu- 
Since B^^-^ is orthogonal to we have 
It follows from Eqs. m and that 

pf^^B’' = = -nxu^Bl^ . ( 21 ) 

Hence we have 



P^^^B’' 

riyU'^ 


( 22 ) 


Evaluating the time and spatial components of Eq. 122) 
gives 


B 


0 

(“) 



UiB'- ja , 
PVa + 


(23) 

(24) 


In the MHD limit, Tg'))) can be expressed simply in terms 
of 6^ as [cf. Eq. 11511 ] 

- b^b'' , (32) 

and the total stress tensor is given by 

= {poh + b'^)u^u'' + (^ + y ) 5 '"'" “ 

The evolution equations for the fluid are given by the 
baryon number conservation equation \/^(pou'^) = 0 and 
the energy-momentum conservation equation = 

0. Conservation of baryon number gives 


dtp^ -I- dj{p^v^) = 0 , (34) 

where p* = p^u^. The spatial components of the 
energy-momentum conservation equation give the mo¬ 
mentum equation 


The evolution equation for the magnetic field can be 
obtained in conservative form by taking the dual of 
Maxwell’s equation = 0. One finds 

V.P*'^" = ^5,(y^P*^") = 0, (25) 

y~9 

where \/—g = . Note that = P*/q^ [see 

Eq. 113ll 1. The time component of Eq. I25II gives the 
no-monopole constraint 

djW = 0 , (26) 

where 

& = ^B^ . (27) 

The spatial components of Eq. (E3 give the induction 
equation 


dtS, + d,{a^Th) = , (35) 

where we have defined the momentum density variable 

s^ = VlS^ = a^T\ 

= {p*h + au°y/yb'^)ui-. (36) 

The time component of the energy-momentum con¬ 
servation equation gives the energy equation. Following 
Font et al [43 , we use the energy variable 

f = - p* = - p* (37) 

The energy equation is then given by 

dtT + - p^v'-) = s , (38) 

where the source term s is 


dtB^ + 9,[V^(w^B(„) - u^Bl^^)] = 0 . (28) 

It follows from Eq. (1^ that 

u^Bl^ - u^Bl^ = {v^B^ - FB^)la , (29) 

where z)® = u^/u^. Hence the induction equation can be 
written as 

dtB^ + dj{v^B^ -F&)=Q . (30) 

C. Evolution of the hydrodynamics fields 

In the literature, a magnetic 4-vector b^ is often intro¬ 
duced. It is related to by 

B^ 

, (31) 

v47r 


s = 

= [(r°°/3*/3^' + 2r°*/3^' -f 

-(r°°/3*-f r°05*a] . (39) 

To complete the system of equations, it remains only 
to specify the equation of state (EOS) of the fluid. In 
this paper, we adopt a T-law EOS 

P = (T - I)poe, (40) 

where T is a constant. We choose the T-law EOS because 
it simplifies some of the calculations, it is applicable to 
many cases of interest, and it is a standard choice for 
demonstrating new computational techniques in the nu¬ 
merical relativity literature. We note that all equations 
in this section apply for any equation of state. General¬ 
ization to a more realistic EOS is not difficult, and we 
plan to use more realistic EOSs in some of our future 
work. 
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To summarize, the evolution equations for the magne¬ 
tohydrodynamic variables are 



III. IMPLEMENTATION 

We use a cell-centered Cartesian grid in our three- 
dimensional simulations. Sometimes, symmetries can be 
invoked to reduce the integration domain. For octant 
symmetric systems, we evolve only the upper octant; for 
equatorially symmetric systems, we evolve only the upper 
half-plane. For axisymmetric systems, we evolve only the 
X — z plane [a (2-|-l)D problem]. In axisymmetric evolu¬ 
tions, we adopt the Cartoon method for evolving the 
BSSN equations, and use a cylindrical grid for evolving 
the induction and MHD equations. 

Our technique for evolving the metric fields is de¬ 
scribed in our earlier papers i mil , so we focus 
here on our MHD algorithms. The goal of this part of 
the numerical evolution is to determine the fundamental 
MHD variables P = {po, , B^), called the “primi¬ 
tive” variables, at future times, given initial values of P. 
The evolution equations iBTIi-iroii are written in conser¬ 
vative form, i.e. they give the time derivatives of the 
“conserved” variables U(P) = (p*, f, 5^, H*) in terms of 
source variables S(P) and the divergence oi flux variables 
F(P): 

atU + V • F = S , (45) 

where F(P) and S(P) are not explicit functions of deriva¬ 
tives of the primitive variables, although they are explicit 
functions of the metric and its derivatives. 

There are several ways of evolving this system. Con¬ 
servative schemes evolve U with the equations gSJ. The 
advantage of this is that highly accurate shock-capturing 
methods can be applied to this set of equations. The 
disadvantage is that, after each timestep, one must re¬ 
cover P by numerically solving the system of equations 
U = U(P), which can be complicated and computa¬ 
tionally expensive. Non-conservative schemes, on the 
other hand, evolve variables which are more simply re¬ 
lated to P but whose evolution equations are not of the 


form of Eq. In such schemes, high-resolution shock¬ 
capturing methods cannot be used and artificial viscosity 
must be introduced for handling discontinuities, but the 
recovery of P is fairly straightforward. After implement¬ 
ing both types of schemes, we found our conservative 
scheme to be more stable when strong magnetic fields 
are present. We evolve Eq. (gni) using a three-step iter¬ 
ated Crank-Nicholson scheme as in several of our previ¬ 
ous papers . This scheme is second order in time and 
will be stable if At < min(Aa:®)/cinax) where in our case 
Cmax is the speed of light. In most applications below, we 
set At = 0.5min(Aa;®). For each Crank-Nicholson sub¬ 
step, we first update the gravitational field variables (the 
BSSN variables). We then update the electromagnetic 
fields by integrating the induction equation. Next, the 
remaining MHD variables (p*, f, and Si) are updated. 
Finally, we use these updated values to reconstruct the 
primitive variables on the new timestep. 


A. The reconstruction step 

We have implemented an approximate Riemann solver 
to handle the flux term in Eq. Below, we will 

demonstrate how the flux / is calculated for a given con¬ 
served variable, u. For simplicity, we will consider the 
one-dimensional case. The generalization to three dimen¬ 
sions is straightforward. The first step in calculating this 
flux is to compute Pl = Pj+i/ 2 _e and P/i = Pi+i/ 2 +e, 
i.e. the primitive variables to the left and right of the grid 
cell interface. We have implemented several methods for 
computing P^ and P/j. 

1 ) Monotonized central (MC) reconstruction 

This method gives second-order accurate results at 
most points. For a given primitive variable p, one sets 

PL = Pi -\-'VpiAxj2 

PR = pi+i - Vpi+iAx/2 . (46) 

Here, Vp is the slope-limited gradient of p: 
Vp = Aa;“^MC(pi+i - p^,pi - Pt-i), where 

MPr W 0 if a 6 < 0 , 

' ( sign(a) min( 2 |a|, 2 | 6 |, |a-I- 6 |/ 2 ) otherwise. 

(47) 

This scheme becomes first-order accurate at extrema of 
P- 

2) Convex essentially non-oscillatory (CENO) recon¬ 
struction 

In this scheme one uses polynomial (usually 

quadratic) interpolation to find cell face values. For 
smooth monotonic functions, these values are accurate 
to third order in Ax. As in the above method, pR and 
Pr usually differ, and the scheme becomes first order at 
extrema ofp. See for details of this reconstruction 

method. 

3) Piecewise Parabolic (PPM) reconstruction 

For smooth monotonic functions on uniform grids, PPM 
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reconstructs face values to third-order accuracy in Ax. 
In this scheme, and pji are equal except in special cir¬ 
cumstances, usually involving shocks M- We have made 
slight modifications to this scheme, the details of which 
are discussed in Appendix^ These changes do not affect 
the fluid evolution in any of the applications below ex¬ 
cept that of unmagnetized stars. In this particular case, 
it is necessary to distinguish between the standard PPM 
scheme as given in |5l) and ours, and so we refer to the 
former as PPM and the latter as PPM+. 

We note that, even when using higher order reconstruc¬ 
tion schemes such as CENO and PPM, our overall evolu¬ 
tion scheme remains second-order accurate in space and 
time. This is because our finite differencing of Eq. gSl) 
is only second-order accurate (although this could be im¬ 
proved), and the BSSN variables are only evolved with 
second-order accuracy. Nevertheless, higher order recon¬ 
struction schemes can provide more accurate results for 
some applications. 


B. The Riemann solver step 


Next, we take the reconstructed data as initial data for 
a piecewise constant Riemann problem, with P = on 
the left of the interface, and P = P/j on the right of the 
interface. The net flux at the cell interface is given by 
the solution to this Riemann problem. 

We use the HLL (Harten, Lax, and van Leer) approx¬ 
imate Riemann solver [H2- The HLL solver is one of 
the simplest shock-capturing schemes as it does not re¬ 
quire knowledge of the eigenvectors of the system. Nev¬ 
ertheless, when coupled to a higher order reconstruction 
method such as PPM, even simpler Riemann solvers, let 
alone HLL, have been shown to perform with an accu¬ 
racy comparable to more sophisticated solvers in shock 
tube problems |49l 1 5, 'll and in binary neutron star simu¬ 
lations [^. To compute the HLL fluxes, one only needs 
to provide a maximum left-going wave speed c+ and a 
maximum right-going wave speed c_ on both sides of 
the interface. Defining Cmax = niax(0, c+a;, c+l) and 
Cmin = — min(0, C-h^C-l), the HLL flux is given by 


j. Cniin/*i? T CjnaxJ^L CminCmax(^i? '^l) f ao\ 

h+1/2 ---——- • (48) 

*-max “T ^min 


We compute the wave speeds c± as described in Section 
3.2 of |3r|. Only the maximum wave speeds in either 
direction along the three coordinate axes are required. 
To determine the speeds in the x direction, one solves 
the dispersion relation for MHD waves with wave vectors 
of the form 


factor of < 2 (thus making the evolution stabler, but also 
more diffusive). In the frame comoving with the fluid, the 
approximate dispersion relation for MHD waves is 

= [va + c^l - v\)\ kl^ , (50) 


where Cg = -y/PP/ (hpo) is the sound speed, va is the 
Alfven speed, and the subscript “cm” refers to comoving 
frame values. To solve this equation for oj/ki in the grid 
frame, we use tUcm = —k^u^, v\ = and k^^ = 

where £ = p^h + 5^ and + u^Uv)k'' is 

the part of the wave vector normal to . Then, for k^, 
one substitutes Eq. lO or its y-axis or z-axis equivalent. 


C. Recovery of Primitive Variables 


Having computed U at the new timestep, we must use 
these values to recover P, the primitive variables on the 
new time level. This is not trivial because, although the 
relations U(P) are analytic, the inverse relations P(U) 
are not. In general, one can do the inversion by numeri¬ 
cally solving a system of nonlinear algebric equations |5(i| . 
Here, we discuss the case where the EOS is given by 
Eq. (Iinil . We need to solve the following four equations 
for the variables e and Ui [c.f. Eqs. (El, (IS3]: 

0 = Pi^hui + aA/^u^b^Ui-ay^b°bi-Si (51) 

0 = {au^ — 1+ reau^)pi, + Ayjb^{au^)^ 

(^+y) (52) 

In this system, the variables h, it°, 6^, 6°, bi, and 
P are treated as functions of the unknown variables e 
and Ui, the known set of conserved variables U, and 
the known metric quantities. The primitive variables 
P = {po,P,v'^,B^) are then constructed using the e and 
Ui which solve the above system. The updated values of 
P® are already known from the induction step. To ob¬ 
tain the remaining primitive variables, the following set 
of steps may be used: 


= 

i(l -1- 

(53) 

Po = 

P* 

(54) 

a^Pyu^ 

P = 

(T - l)poe 

(55) 

P = 

- / 3 ‘ . 

(56) 


k^ = {-co,ki,0,0). (49) 

The wave speed is simply the phase speed oj/ki. The 
speeds along y and z are computed in a similar way. As 
in |3r)|| , we replace the full dispersion relation by a simpler 
expression which overestimates the maximum speeds by a 


As the primitive variable inversion is much simpler 
without magnetic fields, we use a different scheme in such 
cases. First, the condition u^u^ = — 1 is rewritten as 

= pI + , (57) 
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where w = avPPi,. Using the definition of t in Eq. ea, 
one may write h in terms of w and the conserved vari¬ 
ables: 


rw(f- (r- i)p; 
rri ;2 - (r - l)p 2 


(58) 


Substituting this expression into Eq. leads to a quar- 
tic equation for {w — p*) We solve this equation 

using a standard polynomial root finder, and then find h 
by substituting w back into Eq. The primitive vari¬ 
ables can then be constructed according to the following 
set of steps: 


u° = 

w 

(59) 


ap* 

Po = 

pI 

(60) 


Vjw 


r -1 


p = 

-j-po{h-l) 

(61) 

p = 

1 ^ 

_ gi 

Pi,h 

(62) 


D. Low-Density Regions 


If a pure vacuum were to exist anywhere in our com¬ 
putational domain, the MHD approximation would not 
apply in this region, and we would have to solve the 
vacuum Maxwell equations there (see |57| for an exam¬ 
ple). In many astrophysical scenarios, however, a suffi¬ 
ciently dense, ionized plasma will exist outside the stars 
or disks, whereby MHD will remain valid in its force- 
free limit. For the code tests involving magnetic fields 
which we will be presenting in this paper, there is no 
such low density region and no special treatment is re¬ 
quired. We do, however, present tests below with un¬ 
magnetized, rotating stars. For these tests, we do not 
impose floors on the hydrodynamic variables. This is the 
“no-atmosphere” approach used in M- However, in the 
low-density regions near the surface of the star, we some¬ 
times encounter problems when recovering the primitive 
variables; in particular, the equations U = U(P) have no 
physical solution. Usually, unphysical U are those corre¬ 
sponding to negative pressure. At these points, we apply 
a fix, first suggested by Font et al ^5. In the system 
of equations to be solved, we replace the energy equa¬ 
tion with the adiabatic relation P = kpq, where k 
is set equal to its initial value. This substitution gau- 
rantees a positive pressure. When magnetic fields are 
present, the no-atmosphere approach is not suitable, and 
a very small positive density must be maintained outside 
the stars. Special techniques for dealing with the low- 
density region in MHD calculations have been explored 
in 


E. Constrained Transport 

Unphysical behavior may be expected if the divergence 
of the magnetic field is not forced to remain zero. Thus, 
constrained transport schemes have been designed to 
evolve the induction equation while maintaining diB'^ = 0 
to roundoff precision [^. We use the flux-interpolated 
constrained transport (flux-CT) scheme introduced by 
Toth 1^ and used by Gammie et al |33- This scheme 
involves replacing the induction equation flux computed 
at each point with linear combinations of the fluxes com¬ 
puted at that point and neighboring points. The combi¬ 
nation assures both that second-order accuracy is main¬ 
tained, and diB'- = 0 is strictly enforced. 


F. Black Hole Excision 

Black hole spacetimes are evolved using singularity ex¬ 
cision. This technique involves removing from the grid a 
region (the “excision zone”) containing the spacetime sin¬ 
gularity. Rather than evolving inside this region, bound¬ 
ary conditions are placed on the fields immediately out¬ 
side the excision zone. If the region excised is inside the 
event horizon, the causal properties of the spacetime will 
prevent the effects of excision from contaminating the 
evolution outside the black hole. Our excision zones are 
spherical and are placed well inside the apparent hori¬ 
zons, hence well inside the event horizons. For details on 
the excision boundary conditions placed on the metric 
fields, see |^. In we set the hydrodynamic vari¬ 
ables equal to zero at the excision boundary (i.e. matter 
is destroyed when it hits the excision zone). With our 
new code, we find that this excision boundary condition 
for the fluid variables is still adequate in the absence of 
magnetic fields. When magnetic fields are present, how¬ 
ever, it can become problematic. Therefore, we now set 
the MHD variables on the excision boundary by linearly 
extrapolating the primitive variables along the normal to 
the excision surface. 


IV. CODE TESTS 
A. Unmagnetized Relativistic Stars 

In this section, we test the ability of our GRMHD 
code to handle rotating relativistic stars without mag¬ 
netic fields. For initial data, we take a perfect fluid 
with a polytropic equation of state P = with 

n = 1, and we choose our units such that k = 1. (For 
a description of the code used to generate these rotat¬ 
ing equilibrium stars, see [fij- Eqs. (15)-(23) of [fij 
give the scaling relations to arbitrary k.) We evolve 
the three rotating polytropes described in Table D We 
adopt equatorial symmetry in all cases. We note that 
stars A, B, and C in Table ^ correspond, respectively, to 





TABLE I: Rotating Equilibrium Stars {n — 1, k — 1)“. 


Star 

M " 

p c 
-^eq 

p d 

iXc 

7^ " 


r/|wK 

Hp/Heq® 

A 

0.170 

0.540 

0.881 

0.88 

0.35 

0.032 

1.00 

B 

0.171 

0.697 

0.780 

0.87 

0.34 

0.031 

1.00 

C 

0.279 

1.251 

1.613 

0.30 

1.02 

0.230 

2.44 


“ The maximum ADM mass for a nonrotating n = 1, k = 1 
polytrope is Mmax = 0.164. 

ADM mass 

^ coordinate equatorial radius 
circumferential radius at the equator 
® ratio of polar to equatorial coordinate radius 
^ ratio of rotational kinetic to gravitational potential energy 
® ratio of polar (central) to equatorial angular velocity 

stars C, D, and E in Table I of ^3- First, we demon¬ 
strate convergence for axisymmetric evolutions of star 
A, a uniformly rotating, stable star. This star is known 
to be secularly stable by the turning-point theorem [^ . 
and it is known to be dynamically stable from previous 
numerical simulations j4ll |. Thus, we should find that 
the system maintains equilibrium when evolved in our 
code. In Fig.^ we show the error in the central density 
for three short runs with star A at different resolutions. 
This demonstrates that our standard method for hydro¬ 
dynamics (HLL fluxes with PPM'*' reconstruction) leads 
to second-order convergence. 

Evolutions of star A using several different reconstruc¬ 
tion methods are compared in Fig. [3 Reconstruction 
with the MC limiter leads to a downward drift in the 
central rest-mass density (~ 10% in lOProt, where Prot 
is the rotation period). A similar drift has been seen in 
simulations using other codes [fi^ 1^ . We find that the 
drift converges to zero faster than second order as the res¬ 
olution is increased. CENO reconstruction gives a slight 
improvement, but introduces high frequency oscillations 
in the central density. In Fig. |21 these oscillations are 
not individually distinguishable, but instead make the 
CENO line appear thicker than the others. These oscil¬ 
lations are an artifact of the coordinate singularity near 
the axis in cylindrical coordinates. (We do not see the 
oscillations in 3D runs with CENO.) The high frequency 
oscillations can be removed by adding high-order dissi¬ 
pation. We note, however, that the other reconstruction 
methods represented in Fig. 12 do not display high fre¬ 
quency oscillations and thus do not require such fixes. 
The best results are achieved with PPM and PPM+ re¬ 
construction. Much of the drift in the central density 
vanishes when standard PPM is used, which is consistent 
with the result reported in |^. However, with PPM+, 
the drift is eliminated almost entirely. 

Next we check the ability of our code to distinguish 
radially stable from radially unstable stars. We con¬ 
sider two uniformly rotating stars, stars A and B, which 
are members of a constant angular momentum sequence, 
J = 0.01 in our G = c = k = 1 units. The J = 0.01 
sequence has a turning-point at central rest-mass density 



FIG. 1: Relative error in the central rest-mass density for uni¬ 
formly rotating star A. The error is plotted for three axisym¬ 
metric runs, adopting equatorial symmetry (with resolutions 
32^, 64^, and 128^), and the curves are scaled for second- 
order convergence. All runs are performed with the standard 
hydrodynamic scheme (HLL with PPM"*"). Outer boundaries 
are placed at 7.1M. 



FIG. 2: Normalized error in the central rest-mass density for 
star A with different reconstruction methods. All runs are 
axisymmetric and equatorially symmetric with a resolution 
of 64^ and outer boundaries at 7.1M. With PPM, the drift 
in the central density is strongly reduced with respect to the 
MC and CENO results. With PPM*^, the central density drift 
disappears almost entirely. 
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FIG. 3: Axisymmetric evolution of uniformly rotating stars. 
Star A (solid lines) is stable, while star B (dashed lines) is 
unstable to collapse. The upper window shows the central 
density normalized to its initial value, while the lower gives 
the central lapse. The solid dot indicates the first appearance 
of an apparent horizon during the collapse of star B. 


^crit _ Q_ 32 ^ which has the maximum mass Mmax = 0.172 
for the sequence. For a sequence of uniformly rotat¬ 
ing stars, this turning point marks the onset of secular, 
not dynamical, radial instability |62|. but prior numerical 
simulations |4l| have found the point of onset of dynam¬ 
ical instability to be very close to the point of onset of 
secular instability. We pick two similar stars on either 
side of the onset of secular instability: star A with ini¬ 
tial central rest-mass density Pc(0) = 0.24 on the stable 
branch and Star B with Pc(0) = 0.37 on the unstable 
branch. In Fig. |3 we see that the code correctly finds 
star A to be stable and star B to be unstable. 


Star B collapses to a black hole. Without excision, the 
extreme density and spacetime curvature at the center 
of the collapsing star cause the code to crash shortly af¬ 
ter the formation of an apparent horizon which envelops 
some, but not all, of the star. The evolution can be con¬ 
tinued by excising a region inside the horizon. When we 
do this, we find that all of the matter falls into the hole 
within a few M of the time excision is introduced, leaving 
a vacuum Kerr black hole with roughly the same M and 
J as the initial star B. We then continue to evolve for an¬ 
other 30M. We find that the hole’s angular momentum, 
computed as the sum of surface and volume integrals, 
decreases with time, and this angular momentum drift 
limits the length of time that our evolution remains ac¬ 
curate. C omp arable angular momentum loss was also 
present in . Since this drift appears after most of the 
matter has fallen into the excision region, the source of 



0 1 2 

X/M 


FIG. 4: Snapshots of the rest density contours and the veloc¬ 
ity field iy^ ,v^) in the meridional plane during the axisym¬ 
metric collapse of uniformly rotating star B to a Kerr black 
hole. The contour lines are drawn for po = ]^Q-(o.3j+o.09)^Maa: 
for j = 0, 1,.., 12, where p^°‘^ is the maximum of po at the 
time of excision. The thick dashed curve marks the excision 
zone. The thick solid curve is the apparent horizon. We show 
the system at the time of excison (upper panel), at the time 
at which the last of the matter falls into the excision region 
(middle panel), and at a late time just prior to the termina¬ 
tion of the integration (lower panel). 

the error resides in the evolution of the BSSN variables. 
Since we use the same algorithm to evolve the metric 
as in if is not surprising that the J-drift has the 
same magnitude. We simulated the collapse of star B 
on both 64^ and 128^ grids, and we found, as in 
that the J-drift converges to zero with increasing resolu¬ 
tion. In Fig.^] we show sna psh ots from our post-excision 
evolution on the 128^ grid |^. To check that the final 
geometry corresponds to a stationary Kerr spacetime, we 
confirm that the area and the equatorial and polar cir¬ 
cumferences of the apparent horizon agree with the ap¬ 
propriate values for a Kerr black hole with the M and J 
of our spacetime to better than 1% (see for details). 

We now demonstrate the ability of our code to handle 
differential rotation in both 2D and 3D by considering 
the evolution of star C, a stable, hypermassive, differen- 
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FIG. 5: Snapshots of the angular velocity prohle for an evolu¬ 
tion of the differentially rotating star C in axisymmetry (48^) 
with outer boundaries at 7.1M. The prohle is well-maintained 
in the bulk of the star for over 15 central rotation periods 
(Prct). 


tially rotating star. We evolved this star in axisymmetry 
with PPM+ using a fairly low resolution of 48^ zones and 
outer boundaries at 7.1M. Figure El shows the angular 
velocity profile at several times during this evolution, and 
demonstrates that our code correctly maintains the dif¬ 
ferential rotation. Slight errors in the angular velocity 
arise near the origin. The central angular velocity is par¬ 
ticularly susceptible to error, as its calculation involves 
dividing the local azimuthal velocity by the (small) ra¬ 
dius. During the first 15Prot of this evolution (where Prot 
refers to the central rotation period), the central density 
remains within 5% of its initial value and the ADM mass 
is conserved to within 0.5%, even at this low resolution. 
(The angular momentum computed as a volume integral 
over all space is conserved exactly in our code in axisym¬ 
metry.) During this period, the normalized Hamiltonian 
and momentum constraints (see m for definition) are 
~ 1 %. 

Next, we compare the axisymmetric run of star C de¬ 
scribed above with an equivalent 3D run in equatorial 
symmetry. The 3D run was performed with a 96^ x 48 
grid in {x, y, z) and outer boundaries at 7.1M, so that the 
grid cell size is the same as in the 2D run. The deviations 
in the central density for these two runs are compared in 
Fig. El which shows that the axisymmetric and full 3D 
runs have comparable errors. For t < 6Prot, the two runs 
are similar and the angular velocity profiles agree very 
well. After this time, however, the star in the 3D run be¬ 
gins to move away from the center of the grid, eventually 
making contact with the outer grid boundary. This is due 
to accumulated error in the linear momentum and is a 



FIG. 6: Fractional error in the central density versus time for 
evolutions of star C in 2D and 3D. The 3D run was performed 
in equatorial symmetry with resolution 96^ x 48 with outer 
boundaries at 7.1M. The axisymmetric run has the same 
spatial resolution (grid cell size) as the 3D case, and thus 
employs a 48^ grid with boundaries at 7.1M. The errors in 
the central density are of the same order for both runs. 

well known problem associated with evolutions of stars in 
equatorial symmetry |fi(T|| . In our simulations, the effect 
can be reduced by improving spatial resolution. By com¬ 
paring runs of star C on 64^ x 32, 96^ x 48, and 168^ x 84 
grids, we find that the movement of the center of mass 
converges to zero at third-order in spatial resolution. We 
note that the drift can be removed at any resolution by 
employing rr-symmetry, so that the symmetry boundary 
conditions tie the star to the center of the grid. 


B. Minkowski Spacetime MHD Tests 

Komissarov has proposed a suite of challeng¬ 
ing one-dimensional tests of nonlinear, relativistic MHD 
waves in Minkowski spacetime. Most of the tests (except 
the nonlinear Alfven wave test) start with discontinuous 
initial data at x = 0 (see Table , with homogeneous 
profiles on either side. We integrate the MHD equations 
from t = 0 to t = tfinai, where ffinal is specified in Ta¬ 
ble HD for each case. The gas satisfies a F-law EOS with 
F = 4/3. In all the cases, our computational domain is 
X S (—2,2). Our standard resolution is Ax = 0.01 (400 
grid points). We are able to integrate all the cases using 
the MC resconstruction scheme and with a Courant fac¬ 
tor of 0.5. Thus, the number of timesteps for a given test 
is tfinai/(0.5Ax) = 200tfinai- With PPM resconstruction 
scheme, we need to lower the Courant factor to 0.4 for 
the fast shock test. We obtain slightly better results with 
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X 


FIG. 7: Density profiles for the nonlinear wave tests at t = 
tflnai (see Table EJ. Symbols denote data from numerical 
simulations with resolution Aa; = 0.01. Solid lines in the 
upper 6 panels denote the exact solutions |^ . The solid line 
in the last panel denotes a numerical simulation with higher 
resolution, Aa; = 10“^. 


PPM resconstruction scheme for the shock tube tests 1 
and 2. The two resconstruction schemes give comparable 
results for the other tests. Here we present the simula¬ 
tions using the MC resconstruction scheme. Figures EHH 
compare our simulation results (symbols) with the ex¬ 
pected results (solid lines) [b^ . Our numerical results are 
similar to those reported recently for other codes |^l^ . 
Below, we briefly discuss each of the cases we studied. 

Fast and slow shocks In these tests, the initial MHD 
variables on the left [x < 0) and right {x > 0) satisfy the 
special relativistic Rankine-Hugoniot jump conditions for 
MHD shocks As a result, the discontinuity simply 
travels with a certain speed /r without changing its pat¬ 
tern. The fast shock is the most relativistic case of all 
the tests. In the shock frame, the Lorentz factor of the 
upstream flow is vP « 25. The shock moves with a speed 
/i = 0.2. The slow shock is not as strong and it moves 
with a faster speed (/r = 0.5). Our simulations for these 
two tests agree quite well with the exact solutions. In 
the slow shock test, we see small oscillations (due to nu¬ 
merical artifacts) in density on the right side of the shock 
(see Fig. [T)). This numerical artifact is also seen in sim¬ 
ulations with other codes We have performed 

simulations on these two cases using different resolutions 
and found that the errors converge to first order in Ax, 



X 


FIG. 8: Velocity profiles (u^) for the nonlinear wave tests at 
t = tflnai. Symbols denote data from numerical simulations 
with resolution Ax = 0.01. Solid lines in the upper 6 panels 
denote the exact solutions |^ . The solid line in the last 
panel denotes a numerical simulation with higher resolution. 
Ax = 10"T 


which is expected for problems with discontinuities in the 
computational domain. 

Switch-on/ojf rarefaction In these tests, the left and 
right states are connected by a rarefaction wave at t > 0. 
The tests become more challenging when the tangential 
component of the magnetic field (i.e., B^) is switched 
on/off when going from the right state to the left state. 
The exact solutions are obtained by integrating a sys- 
tem of ordinary differential equations (see e.g., l7(l| i. 

Our simulation results agree with the exact solutions very 
well, except that we see numerical artifacts near the trail¬ 
ing edge of the rarefaction wave in the switch-off test. 
We also observe a small oscillation (not visible on the 
scale shown in Figs. 0 and 0) near the leading edge of 
the rarefaction wave in the switch-on test. These nu¬ 
merical artifacts are also seen in simulations with other 
codes [ 33,113 • As explained in |b3| , the oscillation results 
from perturbations created by numerical dissipation dur¬ 
ing the initial stage when the wavefront is very steep. The 
perturbations propagate across the main wave and then 
separate from it. 

Shock tubes 1 and 2 The initial left and right states 
are given in Table im At time t > 0, the left and right 
states are connected by a rarefaction wave, a contact dis¬ 
continuity and a shock wave. The exact solution can be 
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TABLE II: Initial states for one-dimensional MHD tests. 


Test 

Left state 

Right State 

^final 

Fast Shock 
{g = 0.2^=) 

N = (25.0,0.0,0.0) 
B^/y/W = (20.0, 25.02,0.0) 
P= 1.0, po = 1.0 

N = (1.091,0.3923,0.00) 
P77I)/= (20.0,49.0,0.0) 
P = 367.5, po = 25.48 

2.5 

Slow Shock 
{g = 0.5^=) 

u* = (1.53,0.0,0.0) 
B^/y/W = (10.0,18.28,0.0) 
P = 10.0, po = 1.0 

N = (0.9571,-0.6822,0.00) 
P7v^= (10.0,14.49,0.0) 
P = 55.36, po = 3.323 

2.0 

Switch-off Fast 
Rarefaction 

N = (-2.0, 0.0, 0.0) 
B^/y/W = (2.0, 0.0, 0.0) 

P= 1.0, po = 0.1 

u' = (-0.212,-0.590,0.0) 
P7\^= (2.0,4.71,0.0) 

P = 10.0, po = 0.562 

1.0 

Switch-on Slow 
Rarefaction 

N = (-0.765,-1.386,0.0) 
P7v^= (1.0,1.022,0.0) 
P = 0.1, po = 1.78 X 10"® 

u' = (0.0, 0.0, 0.0) 
B^/^/W = (1.0, 0.0, 0.0) 

P = 1.0, po = 0.01 

2.0 

Shock Tube 1 

N = (0.0, 0.0, 0.0) 
B^/^/W = (1.0, 0.0, 0.0) 

P = 1000.0, po = 1.0 

N = (0.0, 0.0, 0.0) 
p774)F= (1.0,0.0,0.0) 
P= 1.0, po = 0.1 

1.0 

Shock Tube 2 

N = (0.0, 0.0, 0.0) 
B'/y/W ^ (0.0,20.0,0.0) 

P = 30.0, po = 1.0 

= (0.0, 0.0, 0.0) 
P7V4^= (0.0, 0.0, 0.0) 
P= 1.0, po = 0.1 

1.0 

Collision 

N = (5.0, 0.0, 0.0) 
B^/y/W = (10.0,10.0,0.0) 
P= 1.0, po = 1.0 

N = (-5.0, 0.0, 0.0) 
P7V4^= (10.0,-10.0,0.0) 
P= 1.0, po = 1.0 

1.22 

Nonlinear Alfven wave“ 
{g = 0.626*’) 

N = (0.0, 0.0, 0.0) 

= (3.0, 3.0, 0.0) 
P= 1.0, po = 1.0 

N = (3.70,5.76,0.00) 
B'/y/W = (3.0, -6.857, 0.0) 
P= 1.0, po = 1.0 

2.0 


In all cases, the gas satisfies the F-law EOS with F = 4/3. For the first 7 tests, the 
left state refers to x < 0 and the right state, a; > 0 . 

^ ft is the speed at which the wave travels 
For the nonlinear Alfven wave, the left and right states are joined by a continnous 
function. See 173 or Appendix IHlfor details. 


computed using a method similar to El In the shock 
tube 1 test, the solution consists of a thin layer of shocked 
gas, which is poorly resolved in our simulation and has 
a wrong value of shell density. The thin layer is covered 
by only 5 grid points with our resolution. We found that 
the correct density is obtained in higher resolution simu¬ 
lations in which Aa: < 0.0035, which provides > 12 grid 
points across the thin layer. Our results in these two 
tests are comparable to those reported in . 

Collision In this test, the flows on both sides travel 
with equal speed but in opposite directions. The tan¬ 
gential component of the magnetic field is also equal in 
magnitude but opposite in direction. We do not have the 
exact solution for this test. Thus, we compare our lower- 
resolution (Aa; = 0.01) simulation with a high-resolution 
one (Aa; = 10“^). The lower-resolution simulation re¬ 
sults are qualitatively the same as the results reported 
in 1^ , but not as good as the results of Komissarov , 
who uses a more sophisticated Riemann solver. 

Nonlinear Alfven Wave The initial data for this test are 
qualitatively different from the other seven tests. The left 
(a; < —W/2) and right (a; > W/2) states are separated by 
a width W = 0.5 at t = 0. The two states are joined by 
continuous functions in the region x G {—W/2, W/2) at 


t = 0. The details of the setup of initial data can be found 
in [t^. which we summarize in Appendix f5l The pattern 
should simply move with a constant speed = 0.626. 
FigureElshows the simulation results (symbols) and exact 
solution (solid lines). The simulation results are again 
similar to |67| . Since there are no discontinuities in this 
problem, we expect the errors to converge at second order 
in Ax. To demonstrate this, we consider a grid function 
g with error 6g = g — We calculate the LI norm 

of Sg (the “average” of 6g) by summing over every grid 
point i\ 

N 

Ll{6g) = AxJ2\9^-9^^^^\^^)\, (63) 

i=l 

where N cx 1/Aa; is the number of grid points. Figure [Tni 
shows the LI norms of the errors in and 

at t = tgnai = 2.0. We find that the errors in 
and B^ converge at second order in Ax. The error in B^ 
converges at slightly better than second order in Ax. 


C. Curved Background Tests: 
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FIG. 9: Nonlinear Alfven wave test. Symbols are simulation 
results with resolution Ax = 0.01 and solid lines are the exact 
solution. The profiles are shown at time t = tfinai = 2.0. 
Our computational domain is a; £ (~2, 2). We only show the 
region 0.3 < x < 2.0 in this graph. 



FIG. 10: LI norms of the errors in u^, and for the 

nonlinear Alfven wave test at t = tfinai = 2.0. This log-log plot 
shows that the LI norms of the errors in and B^ are 

proportional to (Ax)^, and are thus second-order convergent. 
The error in B^ goes as a slightly higher power of Ax. 


Relativistic Bondi Flow 

Next, we test the ability of our code to accurately 
evolve the relativistic MHD equations in the strong grav¬ 
itational field near a black hole. Specihcally, we check 
its ability to maintain stationary, adiabatic, spherically 
symmetric accretion onto a Schwarzschild black hole, in 
accord with the relativistic Bondi accretion solution . 
It has been shown that the relativistic Bondi solution is 
unchanged in the presence of a divergenceless radial mag¬ 
netic field 1^. The advantages of this test are that it 
involves strong gravitational helds and relativistic flows, 
and that there exists an analytic solution with which to 


test our results. We write the metric in Kerr-Schild (in¬ 
going Eddington-Finklestein) coordinates; in this way, all 
the variables are well behaved at the horizon (“horizon 
penetrating”), and the excision radius can be placed in¬ 
side the event horizon. We begin by holding the metric 
field variables hxed in order to prevent the black hole 
from growing due to accretion. With the metric fixed, 
the flow is exactly stationary in the continuous limit. 
When evolved with a finite-difference code, discretization 
errors will cause small deviations in the flow from its ini¬ 
tial state. These deviations should converge to zero as 
resolution is increased. Eventually, the system may set¬ 
tle down to an equilibrium solution of the descretized 
hydrodynamic equations. To diagnose the behavior of 
our code, we introduce two variables. The deviation of 
the fluid conhguration from the analytic Bondi solution 
we measure by the LI norm in 3D of |p* — 

(where is the analytic value of p*), normalized by 

the rest mass: 

^ ^xAyAz yj,Zk)\ 

AxAyAzY..^^j,pf^^^*{xi,yj,Zk) 

To measure the settling down of the solution to a numer¬ 
ical equilibrium, we monitor Ap;,, the L2 norm of pbAt: 


1/2 


Apb = AtAxAyAz 


i,j,k 


(65) 


The quantities dp* and Apt were chosen because they 
correspond to the diagnosics used to monitor Bondi ac¬ 
cretion test problems in and ^3: respectively. 

For this test, we evolve the same configuration 
used by 1?^. I 23 , and [ 33 . The sonic radius is at 
Schwarzschild (areal) radius Vg = 8M, the accretion rate 
is M = 1, and the equation of state is T = 4/3. As in 
we set our excision radius at rex = 1.9M (the horizon is 
at 2M) and evolve for lOOM. (We hnd that the system 
settles to equilibrium long before lOOM.) We place outer 
boundaries at lOM, at which point the analytic values of 
the conserved variables are imposed. 

First, we evolve this accretion flow in the absence of 
a magnetic field. In Figure ^2 we show the results 
for an axisymmetric grid of 64^ using various numeri¬ 
cal techniques. We find that using MC reconstruction 
gives much better results than using PPM. This is prob¬ 
ably due to the larger numerical diffusivity in the MC 
scheme which stabilizes spurious numerical oscillations. 
PPM can be “corrected,” however, by adding a small 
dissipation. In [t^ , the oscillations are removed by shift¬ 
ing the numerical spatial stencil in supersonic flow. We 
have instead addressed the problem by adding a small 
Kreiss-Oliger dissipation |^. We find that PPM with 
Kreiss-Oliger dissipation performs as well as MC for this 
problem (see Fig. EJ. We have also evolved an accretion 
flow with small M while allowing the metric to evolve. 
We find that the flow is stable, and that, as expected, the 
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FIG. 11: Unmagnetized Bondi accretion onto a Schwarzschild 
black hole. Four different methods are compared: MC recon¬ 
struction, CENO reconstruction, and PPM reconstruction, 
both with and without Kreiss-Oliger (KO) dissipation. Each 
run was performed on a 64^ grid using axisymmetry. 


FIG. 12: after lOOM for various initial values of 

b^Polr —2M- PPM with Kreiss-Oliger dissipation is used in 
each run. Results from 64^ and 128^ grid runs are compared. 
The (5p* found with the 128^ grid is multiplied by 4 to allow 
scaling to be checked for second order convergence. 


irreducible mass of the black hole slowly increases [t^ . 
The MHD variables remain near the Bondi equilibrium 
initial values until the black hole grows appreciably. 

Next, we evolve with a radial magnetic field. In 
the continuum limit, the magnetic forces cancel exactly. 
However, the cancellation will not be exact in a finite- 
difference code, and this test can be quite difficult for a 
GRMHD code when the magnetic field is strong. 

In Fig. El we plot the error, measured by 6p^, after 
lOOM of evolution, for 2D (axisymmetric) runs with var¬ 
ious values of b'^/po at the horizon. We use the PPM 
reconstruction method with Kreiss-Oliger dissipation for 
each run. In order to test convergence, we use both 64^ 
and 128^ grids. We are able to evolve with magnetic fields 
b'^ /Po\r= 2 M ^ 30. Stronger radial fields quickly crash 
the code. For &^/po|r= 2 M ^ 5, we find that settles 
quickly to a final value. For larger &^/po|r= 2 M, <5p* does 
not settle as well, but the results are still second-order 
convergent after lOOM. Evolving with MC gives better 
“settling” behavior for 5 < fe^/po|r= 2 M ^ 30, but it is still 
not possible to evolve flows with b"^/ po\r= 2 M ^ 30. Gam- 
mie, McKinney, and Toth are able to evolve with 
much higher 6 ^/po|r= 2 M- This is probably due to their 
use of a spherical-polar coordinate grid, which is better 
adapted to the spherical symmetry of this problem than 
our cylindrical grid. 

We have also evolved the Bondi flow in three dimen¬ 
sions on 64^ grids using octant symmetry. We find that 
we can maintain equilibrium flow for 6 ^/po|r= 2 M ^ 10 in 
3D. 


D. Dynamical Background Tests: Gravitational 
Wave-Induced MHD Waves 

To test the capability of our code to handle dynami¬ 
cal gravitational and MHD fields simultaneously, we con¬ 
sider a gravitational wave oscillating in an initially homo¬ 
geneous, uniformly magnetized fluid. The gravitational 
wave will, in general, induce Alfven and magnetosonic 
waves In a companion paper (hereafter. 

Paper H), we perform a detailed analysis of this problem 
and provide an analytic solution for the perturbations in 
a form which is suitable for comparison with numerical 
results. 

This test problem is one-dimensional. We consider 
a linear, standing gravitational wave whose amplitude 
varies in the z-direction: 

h+(t,z) = /i+o sinkz coskt , ( 66 ) 

hx{t,z) = hxosinfczcosfct , (67) 

where k is the wave number, and h+g and hxo are con¬ 
stants. We assume that at t = 0, the magnetized fluid is 
unperturbed: 

P( 0 , z) = Po , po{0, z) = po , ( 68 ) 

P(0,z)=0, B\0,z) = B^^. (69) 

Subsequently, the gravitational wave excites the MHD 
modes of the fluid. As discussed in Paper H, the grav¬ 
itational wave is unaffected by the fluid to linear order, 
and the metric perturbation, z), in the transverse- 

traceless (TT) gauge can be calculated from Eqs. lIHHll 
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and The perturbations in pressure SP{t,z), veloc¬ 
ity Sv'^{t, z), and magnetic field SB^{t, z) can be computed 
analytically as shown in Paper II. Our analytic solutions 
are valid as long as we are in the linear regime in which 
the following three inequalities hold (see Paper II): 


\h^''\ ^ ho 

W^\ ^ 'ho 


< 1 , 

< , 

~ 1 /v^, 


(70) 

(71) 

(72) 


where /iq = y ^ = Po(l + eo) + 7b + ^o- 

The analytic solution is a superposition of the three 
eigenmodes of the homogeneous system (the Alfven, slow 
magnetosonic, and fast magnetosonic waves) and a par¬ 
ticular solution which oscillates at the frequency of the 
gravitational wave. The induced Alfven wave obeys the 
dispersion relation 


=uj\ = {k- vaY , (73) 


where k = kz is the wave vector associated with the 
standing gravitational wave, and va = BojV 47 r£ is the 
Alfven velocity. This mode gives rise to a velocity per¬ 
turbation 


Sv (X iiA = k X Va- (74) 

The frequencies of the induced slow and fast magne¬ 
tosonic modes, ujmi and are found by solving the 

following dispersion relation for 

uj^ - [k'^c^ + cl{k ■ va)‘^]oj'^ + k‘^cl{k ■ vaY = 0 , (75) 

where = Va + Cg(l — Va)- For the corresponding 
eigenvectors, one has: 

6v oc iimi = VA+ , —7^ 1 = 1,2. (76) 

Note that iiA is orthogonal to iimi and Um 2 j but 
and Um 2 are not, in general, orthogonal to each other. 

The setup for our code test is as follows. Our compu¬ 
tational domain is z S (—1,1). We choose k = 2 tt so that 
our computational domain covers two wavelengths of the 
gravitational wave. Our standard resolution is Az = 0.01 
(200 grid points). At time t = 0, we assign the metric 
gfj.ui0, z) = + hf_,^{0, z), where = diag(-l, 1,1,1) 

is the Minkowski metric, and the nonzero components of 
hf^^{0,z) are: 


hxxi0,z) = -hyy{0,z) = h+{0,z) , (77) 

hxyi0,z) = hyx{0,z) = hx{0,z) . (78) 

We choose geodesic slicing and zero shift (a = 1, /3* = 0) 
as our gauge conditions. Hence we set the initial extrin¬ 
sic curvature to zero (771^ ( 0 , z) = 0 ) in accord with our 
gauge choices and Eqs. inSl), (EJi and ©• We also set 


the MHD variables at t = 0 according to Eqs. and 
ra . We choose the adiabatic index P = 4/3 in all of 
our simulations in this section. Periodic boundary con¬ 
ditions on both matter and gravitational field quantities 
are enforced at the upper and lower boundaries in z. We 
expect that the metric, as well as the MHD quantities, in 
our full GRMHD simulations should agree with the an¬ 
alytic solutions given in Paper H to linear order as long 
as the inequalities in Eqs. (COl-IZl are satisfied. 


1. A General Example 

We first consider a general case in which all three of 
the MHD modes are excited. We take the following initial 
data: 


po = 2.78 X 10"® , Po = 1.29 x 10“® , 

Bi = (1.09,8.26,14.4) x lO”'^ , 

h+o = 7x0 = 1-18 X 10-'‘ . (79) 

Figure El gives a comparison of the analytic and numer¬ 
ical solutions for three selected perturbed variables. The 
perturbations are plotted with respect to time for a cho¬ 
sen location on the grid (z = 1/8). Good agreement 
is shown between the numerical and analytic values for 
many periods of the gravitational wave. We also find very 
good agreement for the metric quantities g^x and gxy in 
our simulation and the analytic values calculated from 
Eqs. lIHHIl and The pressure perturbation, however, 
differs from the analytic solution by a slight secular drift. 
(In fact, all variables eventually exhibit a drift away from 
the analytic solution, but the drift is first noticeable in 
the case of the pressure.) 


This secular drift is not due to numerical error, but 
rather is an effect of the nonlinear terms which are ne¬ 
glected in our analytic solution. To show that it is not a 
numerical error, we performed simulations at resolutions 
of 50, 100, and 200 grid points, and we found conver¬ 
gence to second order to a solution with nonzero drift. 
Since the discrepancy is due to nonlinear terms, choosing 
smaller (larger) initial mass-energy density and smaller 
(larger) gravitational wave strength leads to a smaller 
(larger) discrepancy with the analytic solution. In par¬ 
ticular, the size of the discrepancy is controlled by the 
degree to which the conditions in Eqs. fTn ii -fT ^ are sat¬ 
isfied. By evolving a range of initial data sets in which 
we independently varied ho and £/ho, as well as initial 
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FIG. 13: Analytic and numerical solutions for the perturba¬ 
tions of a magnetized fluid due to the presence of a gravita¬ 
tional wave (see Section Ft VD 111 . The thick solid and thin dot¬ 
ted lines represent, respectively, the analytic and numerical 
solutions, though the two lines are not readily distinguishable 
in plots (a) and (b). All quantities are evaluated at z = 1/8 
and are normalized as indicated. Time is normalized by the 
gravitational wave period. 

data sets with /iq = 0 (no gravity wave), we found that 
the numerical solution for the pressure perturbation (SP) 
is always well fit by the relation 

rp _ / dPana + hoPo{ciho + C 2 SIho)t‘^ {ho ^ 0) 

1 c^PoSt^ {ho = 0 ) , 

(80) 

where t is the coordinate time, dPana is the analytic so¬ 
lution for the pressure perturbation given in Paper II, 
and Cl and C2 are constants (see Fig. [Till. (Note that the 
coefficient hoPo is simply the typical scale of the pres¬ 
sure perturbation.) The term proportional to ci corre¬ 
sponds to nonlinear effects of the gravitational wave on 
the fluid, while the term proportional to C2 is related to 
the self-gravity of the fluid |^. Neither of these effects 
are accounted for in our analytic solution [see Eqs. izni- 
m]- Thus, the disagreement of the numerical and ana¬ 
lytic results may be reduced by choosing initial data with 
smaller ho and £/ho- 

To extract the various MHD modes in this test, we 
have performed two projections of the velocity. Fig¬ 
ure [EJr shows the numerical and analytic values of the 
projection along ua (again evaluated at z = 1/8). Be¬ 
cause we have projected the velocity along the direction 
of the Alfven mode eigenvector, which is orthogonal to 
the fast and slow mode eigenvectors, we only pick up 
contributions from the Alfven wave and the particular 



FIG. 14: Relative error in the pressure perturbation (eval¬ 
uated at 2 = 1/8) for three different initial data sets. The 
solid lines give the numerical results while the dashed lines 
give the fits from Eq. if^ . The three initial data sets 
are derived from the standard case in Eq. ifTHli by tak¬ 
ing {po, PoV^'^ = Hpo, {B{,r- = VC(Po)°'^ and 

{/i+o, hxo}"'** = C{^-l-o, where ^ and ( are con¬ 

stants and “old” refers to the values in Eq. I79II . (Effectively, 
r;®"' = and /ig®"' = C^o''^-) The regime of validity for 

the analytic solution is approached for small ho and |T^i/|/ho 
(see Eqs. EU and ED), or equivalently, for decreasing and 
For the curves shown, these values are: (a) = 3, 

l/C = 3, (b) C = 2, 5/C = 2, and (c) C = 1, C/C = 1- Moving 
from (a) to (c), we find that the relative error decreases as 
expected and that the errors are well fit by Eq. I8t)l . Note 
that the normalization differs in the three cases to reflect the 
differing values of Po and ho- 


solution (see Paper II for the analytic expression of the 
particular solution). The lower panel shows these in¬ 
dividual contributions along with the total. One can 
see that both the Alfven and particular components con¬ 
tribute significantly to the velocity perturbation in this 
direction. Next, Fig. (Tfil shows the projection of the ve¬ 
locity along the direction of the slow mode eigenvector. 
This time, there are contributions from the slow and fast 
modes in addition to the particular solution, and one 
again sees that all modes are contributing strongly. (The 
slow and fast modes are both present in this velocity 
projection since Umi and Um 2 are not orthogonal.) Fig¬ 
ures ESI and ESI taken together, show that our code cor¬ 
rectly manifests all three MHD waves, in addition to the 
particular solution contribution. 
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TABLE III: Evaluation of Methods 


Method 

Equilibrium Stars 

Characterization 

Shocks 

Alfven Waves 

HLL-tPPM/HLL-tCT 
HLL-bCENO/HLL-tCT 
HLL-tMC/HLL-bCT 
HLL-tMC/VP 
non-conservative/CT 

performs well 
central density drifts 
central density drifts 
central density drifts 
performs well 

performs well 
performs well 
performs well 
unacceptable oscillations 
problems for high /po 

performs well 
performs well 
performs well 
acceptable 

problems for high f po 



0 2 4 6 8 


t/{27r/k) 


FIG. 15: Velocity projected along the direction of the Alfven 
mode eigenvector {5va = 5v ■ua/\'ua\), for the case discussed 
in Section flV D ll fthe general case), (a) Analytic and numer¬ 
ical plots of the projected velocity. (The curves cannot be 
distinguished on this scale.) (b) Contributions to the ana¬ 
lytic solution for SvA/ho. The particular solution (dashed) 
and Alfven wave (dotted) contributions added together give 
the total perturbation (solid line). All quantities are evalu¬ 
ated at 2 = 1/8. 


2. Special Cases 

To further test our code, we consider several initial 
data sets with special properties. Starting with the ini¬ 
tial data set given in Eq. (EH), we hold £ and /iq con¬ 
stant but change the balance between the plus and cross 
modes so that the magnetosonic modes are not excited. 
As explained in Paper II [in particular, see Eq. (93) and 
surrounding discussion], this occurs when h+o and ft-xo 
satisfy the equation 

[{Blf - - 2B^,Blh^o = 0 . (81) 



FIG. 16: Velocity projected along the direction of the slow 
magnetosonic mode eigenvector (Svm 2 = 5v ■ 'Um 2 /|iim 2 |), 
for the same case as in Figure EH (a) Analytic and numer¬ 
ical plots of the projected velocity, (b) Contributions to the 
analytic solution for 6 vm 2 lho. The fast magnetosonic (dot¬ 
ted), slow magnetosonic (short dashed), and particular so¬ 
lution (long dashed) contributions added together give the 
total perturbation (solid line). All quantities are evaluated at 
z = 1/8. 

Thus, we obtain the new initial data set: 

po = 2.78 X 10“® , Po = 1.29 x 10“® , 

B^q = (1.09,8.26,14.4) x 10“® , 

h+o = 4.31 X 10“® /ixo = 1-61 X 10-^ . (82) 

In the analytic solution for this case, the pressure per¬ 
turbation vanishes identically. In accord with this, our 
numerical solution for the pressure perturbation shows no 
oscillations, though the slight secular drift is still present. 
Similarly, the projection of the velocity along the slow 
mode eigenvector vanishes up to some small-amplitude 
noise, but the projection along the Alfven mode eigen¬ 
vector does not vanish and the analytic and numerical 
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solutions agree very well. Thus, by changing only the rel¬ 
ative proportion of the plus and cross modes, we have ar¬ 
rived at a very different physical outcome from the more 
general case described in Section irvTm and our code 
again correctly identifies the modes which are present. 

The analytic solutions in Paper II also indicate that 
the gravitational wave has no effect on the fluid if (1) 
Bq = Bq = 0, or (2) Bq = 0 and Bq and Bq satisfy 
Eq. inmi . We have performed simulations in these two 
special cases and found that our numerical solutions for 
the perturbations contain only small amplitude noise or 
the secular drift due to nonlinear effects, as expected. 


V. CONCLUSIONS 

We have developed the first code which is able to evolve 
the full coupled Einstein-Maxwell-MHD equations in 3-1-1 
dimensions without approximation. Our code is able to 
model the behavior of magnetized, perfectly-conducting 
fluids in dynamical spacetimes. We have confirmed the 
ability of this code to accurately simulate unmagnetized 
hydrodynamic stars, MHD shocks, Alfven waves, mag¬ 
netized accretion onto a black hole, and the excitation 
of MHD modes in a magnetized fluid driven by gravita¬ 
tional waves. We have performed 1,2, and 3 dimensional 
tests. 

We have tested several different integration schemes. 
In Table Iml we evaluate the behavior of each method 
under various tests. The first row describes the re¬ 
sults obtained by using our “standard” code (listed 
as HLL-I-PPM/HLL-I-CT), in which the fluid equations 
are evolved using HLL fluxes and PPM reconstruction, 
while the magnetic induction equation is evolved us¬ 
ing HLL fluxes interpolated to preserve the constraints 
(CT). The next two rows describe schemes which differ 
from HLL-I-PPM/HLL-I-CT only in the reconstruction 
method used with the hydrodynamic variables (p*, f, and 
Si): HLL-I-CENO/HLL-I-CT uses CENO reconstruction, 
and HLL-I-MC/HLL-I-CT uses MC reconstruction. We 
have also experimented with more significant changes. 
HLL-I-MC/VP evolves the fluid variables with HLL fluxes 
and MC reconstruction (like HLL-I-MC/HLL-I-CT), but 
the induction equation is solved by evolving the vec¬ 
tor potential HL In this method, the magnetic field 
is automatically divergence-free, because it is calculated 
as the curl of the vector potential: S® = Ak^i- 

Finally, we test a nonconservative MHD scheme, non- 
conservative/CT, which is a straightforward extension 
of our previous hydrodynamics code to MHD. Non- 
conservative/CT uses flux-CT to maintain constrained 
transport. It differs from the non-conservative MHD 
code of 1^ in that our grid is not staggered, our en¬ 
ergy variable is different, and our time integration is 
done differently. From the table, we draw several conclu¬ 
sions. (1) PPM is the best of the reconstruction meth¬ 
ods considered here. Using PPM reconstruction, we are 
able to achieve accurate evolutions, even with a simple 


(HLL) Rieman solver. (2) Evolving the magnetic field 
directly with constrained transport gives better results 
than evolving a vector potential. The vector potential 
method performs especially poorly in the presence of 
shocks. (3) Our nonconservative method works well for 
problems involving only weak magnetic fields, but it be¬ 
comes unstable for problems in which the magnetic en¬ 
ergy density significantly exceeds the gas energy density. 
It is, therefore, not suitable for some problems. (Note, 
however, that better nonconservative MHD codes have 
been developed |29| which can accommodate larger mag¬ 
netic fields.) 

Our MHD code has limitations similar to those of other 
MHD codes in the literature. In particular, accurate 
evolution is difficult when ^ po- This could poten¬ 
tially cause problems in the low-density regions in some 
applications. However, the experience of other numeri¬ 
cal MHD groups suggests that these difficulties are sur¬ 
mountable. 

Having passed the tests described above, we will next 
apply our code to the study of self-gravitating magne¬ 
tized fluids in astrophysical problems. In particular, we 
plan to model the braking of differential rotation de¬ 
scribed in several applications in the introduction. We 
also plan to simulate gravitational collapse in order to 
explore the behavior and influence of magnetic fields on 
supernovae and collapsars. 
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APPENDIX A: IMPLEMENTATION OF PPM 
RECONSTRUCTION 

The piecewise parabolic method (PPM) is an algo¬ 
rithm used to construct the values of a primitive variable 
p to the left and right of each zone interface (p^i /2 “ 

Pi+i/ 2 _e and pfj^i /2 ~ Pi+i/ 2 +t)- It consists of several 
steps. First, one interpolates to Pi+ 1/2 according to 

Pi+ 1/2 — Pf+ 1/2 — Pi+l /2 (A-1) 

= P^ + ^{Pi+1 - Pi) + ^Aa;(Vpi - Vpj+i) , 
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where Vp is the MC slope-limited gradient of p (see 
Eq. 14711 1. Note that the factor i on the last term differs 
from the g sometimes appearing in the literature. We 
find more accurate results with the | in Eq. m- Next, 
Pi+i /2 ^'^^Pi+i /2 adjusted using “steepening”, “flat¬ 
tening”, and “monotonizing” algorithms, which are in¬ 
tended to stabilize the evolution and sharpen shock pro¬ 
files. There are several adjustable parameters in the PPM 
scheme; we use the values recommended in |5l). 

As originally proposed, PPM reconstruction reduces 
to first-order accuracy at extrema of p. This is due to a 
“monotonization” step of the PPM algorithm which re¬ 
moves local extrema in the interpolation function in order 
to suppress unphysical oscillations near shocks. Near a 
maximum or minimum of p, however, the extremum in 
the interpolation function represents the true behavior 
in p. Therefore, we have followed [s^ in distinguish¬ 
ing between local extrema caused by numerical oscilla¬ 
tions and physical extrema. When the first derivative 
Aa;“^(pi+i — Pi) changes sign, but the second derivative 
does not change sign over two grid cells in either direc¬ 
tion, the extremum is regarded as a physical maximum 
or minimum and the standard monotonization routine is 
not applied. Third-order accuracy may still be sacrificed 
at these points by other steps of the PPM algorithm, 
but we have found that this modification significantly 
improves accuracy for the evolution of stars. When in¬ 
terpolating poj we also turn off the monotonization when 
Pq is within 15% of its maximum value on the grid. We 
do this to make sure we do not lose accuracy at the cen¬ 
ters of our stars, at which the density is usually a global 
maximum. 

Because the Pf+ 1/2 equal (and hence 

Ur = ul) at most cell interfaces, PPM usually picks up 
no dissipation from the (ur — ur) term in the HLL flux 
formula (Eq. iBHll l. Usually, this is a good thing. For 
cases where some extra dissipation is desirable, such as 
in the relativistic accretion test, we add a small Kreiss- 
Oliger dissipation. This takes the form 


^ (AA:ArAZ)4/3 

dtu--" Cko IQ/S,T ’ 


(A2) 


where is the flat-space Laplacian. We have found 
good results with Cko ~ 0 . 1 . 


APPENDIX B: SETUP OF INITIAL DATA FOR 
THE NONLINEAR ALFVEN WAVE 


One exact solution of the MHD equations is an Alfven 
wave traveling in Minkowski spacetime. The solution was 
derived by Komissarov |72| . Here we summarize the re¬ 
sults. Suppose that the position of the wavefront is given 
by the phase function 


In the Lorentz frame of interest, let p, be the wave speed 
and n be the unit three-vector in the direction of prop¬ 
agation of the wave front. Then, with the appropriate 
scaling of $, we define 

= (- 7 t,n) . (B2) 

Note that <i)Q is proportional to the 1-form ka = (—w, fc), 
which is dual to the propagation 4-vector k°‘. We define 
the following scalars 

a = 

B = 

£ = poh + b"^ (B3) 

The wave speed p can be computed from the equation 
(see, e.g., Eq. (23) of [b^]) 

£a'^-B^ = 0 . (B4) 

Consider two arbitrary points r 4 . and r_ connected by a 
simple Alfven wave, then (see |7^ 1 

[P] = [p,] = [b^] = [p]=o^ 

[a] = [B] = 0, 

[«“] = ^[bl , (B5) 

where [/] = /(r+) — /(r^). Below, we will write / = 
/(r+) and /_ = /(r_). 

Thus, the Alfven wave perturbation at r+ is specified 
by one 4-vector, [6“] or [u“]. This 4-vector has one freely 
specifiable degree of freedom (corresponding to the am¬ 
plitude), the other three being removed by the following 
constraints: 


U°Ua 

= - 1 , 



u°‘ba 

= 0 , 



[a] 

= 0 . 


(B 6 ) 

Substituting -I- [«“] and 6 “ 

that these can be rewritten 

= b°^ + f [u“] 

, we see 

2 m“K] + 

[u“] [Ua] 

= 0 , 


a 

U“)[ltQ,] 

= 0 , 




= 0 , 

(B7) 


One solution for two states connected by an Alfven wave 
is given by Komissarov |72|| . 

The properties of nonlinear Alfven waves are easiest 
to study in the wave frame {p = 0). We denote the 
quantities in this frame by a prime. If we choose the x- 
axis to be normal to the wave front and assume = 
r]^^, then in this frame a = u'“, B — b'^, and we have 


(Bl) 


[w'1 = 0, [r] = 0 , 


$(a;“) = 0 . 
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X = , 

7“] = Xib'l ■ (B 8 ) 

The divergence condition on B'* requires that [B'“] = 0. 
Solving the other constraints, one finds that the trans¬ 
verse components of lie on the ellipse 

a-iib'y + (ai2 -I- a2i)b'yb'^ -\- a22b\ 


3. Choose Bl_ and on the left side. Without loss 
of generality, one may set Bi = 0 . 

4. Calculate from Eqs. II22II and ISU. 

5. Calculate B, a, £ and the wave speed /r from 
Eqs. IB311 and IIB4II . 

6 . Compute b'^ and u't by boosting b^ and to 
the wave frame. Then compute B'^_ from b'^ and 


+ (ai3 + o,3i)by + (023 -l- 032)^2 + 033 — 0 ,(B9) 


where 


flu — 1 — Qy , 

033 = — (c^ -I- d) , 
fli3 = 031 = —cay , 
and where 

u'y_ - ^b'l_ 

- xb'- ’ 

^ Xb"^ 

^ u'°_ - xb'- ’ 

The center of this ellipse is at 


022 = 1 ~ Oj, , 

012 = 021 = —ayUz , 

023 = 032 = -CO 2 , (BIO) 


02 = 


u' - xb'- 
d = b^-b'l . 


where 


D = 


b'l - b^u'l 

B ^2 


(Bll) 


{blX) = {^){<^v.o..) , (B12) 


(B13) 


It is convenient to rewrite the equation of the ellipse 1B9II 
in terms of a free parameter 9 defined by 

b'^ = bl + byz{9) cos 9 , (B14) 

b'" = bl + byz{9)sm9 . (B15) 

Substituting this into Eq. 1B9II . one obtains 


byM = 


d + c^/D 


ail cos^ 9 -I- 2 ai 2 sin 9 cos 9 + 022 sin"^ 9 


(B16) 

One can construct an Alfven wave (propagating in x- 
direction) connecting a left state and a right state as 
follows: 


1. Choose the width W that connects the left and 
right state. In our test, we choose W = 0.5. 


7. Set u"''{x) = it'“ and b"''{x) = 6 '^ (constant every¬ 
where). 

8 . Compute y, Oy, Oz, c, d, a^. 

9. Compute d)), h^. 

10. Choose 9{x) consistent with {b'^_,b'^_). In this pa¬ 
per, we choose [s^ 


9{x) = < 


9i + A sin^ 
, 9i + A 


■k{x+WI1) 

IW 


X < -WI2 
-W/2 <x< W/2 , 
X > W/2 

(B17) 

where A is a freely specifiable constant (“ampli¬ 
tude” of the wave) and d/ is given by 


9i = tan 


-1 


b'- - bl 

yy - bl 


(B18) 


One can verify that our choice of 9(x) gives the 
correct left state, and the right state is determined 
by the value of A. We choose A = tt in our Alfven 
wave test in order to compare our numerical results 
with those by Komissarov . This means that the 
tangential component of b"' is rotated by tt when 
going from the left state to the right state. 

11. Compute b'y{x) and b'^{x) from Eqs. 

fTTlTt and frrrm . 

12 . Use [m'“] = x[5'“] to set {u'^,u''^) as a function of 

X. 


13. Compute u'^{x) and b'^{x) from the relations 
u'^’u'y, = —1 and b'^'u'y = 0 . 

14. Compute b^{x) and m''(x) by boosting b"'‘{x) and 
it'^(x) back to the original frame. 

15. Compute B*(x) from m^(x) and b^{x). 


2. Choose po and P, which are constant throughout 
the wave. 


We use this recipe to construct the initial data for our 
nonlinear Alfven wave test. 
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